#include <cassert>
#include <iostream>

#include "../src/kl.hh"
#include "../src/logistic.hh"

namespace 
{
  using namespace playerpiano;

  struct LogPTestData
  {
    struct
      {
        float w[17 * 5];
        Feature f[4];
      }                 inputs;
  
    float logp[6];
  };
  
  void
  test_logp (void)
  {
    LogPTestData data[] = 
    {{{{0.233259, 0.439907, 0.795762, 0.521043, 0.855965, 0.562346,  0.887541, 0.40101, 0.1029, 0.336449, 0.907054, 0.51278, 0.713948,  0.973245, 0.19742, 0.777745, 0.210109, 0.661834, 0.816352,  0.100476, 0.379556, 0.507045, 0.83227, 0.834844, 0.146297,  0.0671381, 0.0365076, 0.313801, 0.290333, 0.504792, 0.148966,  0.912791, 0.187433, 0.168343, 0.241913, 0.400012, 0.473485,  0.195098, 0.0444922, 0.622266, 0.263376, 0.533264, 0.22814,  0.521791, 0.88382, 0.0262192, 0.395871, 0.686947, 0.737523,  0.959081, 0.359363, 0.373146, 0.44719, 0.454289, 0.210397,  0.460355, 0.259757, 0.285946, 0.968485, 0.060343, 0.786272,  0.0908482, 0.923992, 0.438077, 0.522895, 0.557584, 0.695852,  0.916286, 0.639075, 0.531365, 0.299981, 0.229339, 0.901553,  0.572284, 0.940618, 0.856193, 0.454362, 0.117995, 0.730221,  0.395839, 0.194606, 0.832049, 0.761736, 0.335496,  0.408334}, {{14, 0.7412}, {10, 0.837744}, {0, 0.897419}, {10,  0.885438}}}, {-2.25867, -2.10183, -1.15664, -1.635, -1.48169,  -3.30959}}, {{{0.183616, 0.141892, 0.981133, 0.246363, 0.652251,  0.841911, 0.751794, 0.34481, 0.0799671, 0.901293, 0.895601,  0.890448, 0.961972, 0.171072, 0.499762, 0.695842, 0.129924,  0.409336, 0.164266, 0.287509, 0.388723, 0.571592, 0.266847,  0.40207, 0.205107, 0.4297, 0.285714, 0.155707, 0.552856, 0.587789,  0.53392, 0.810897, 0.472889, 0.686496, 0.63832, 0.920449,  0.510917, 0.515424, 0.138558, 0.224606, 0.380993, 0.106088,  0.974292, 0.937098, 0.99227, 0.534497, 0.707445, 0.535027,  0.787162, 0.104797, 0.421731, 0.37932, 0.234306, 0.517008,  0.887811, 0.568424, 0.761417, 0.830512, 0.249491, 0.647975,  0.2505, 0.315088, 0.110934, 0.423369, 0.869507, 0.208999,  0.136642, 0.486271, 0.877237, 0.674502, 0.429197, 0.951244,  0.090075, 0.569706, 0.00746623, 0.571923, 0.855769, 0.0526975,  0.119655, 0.00349966, 0.0943518, 0.222186, 0.870164, 0.355525,  0.843852}, {{15, 0.907098}, {5, 0.759231}, {6, 0.932156}, {10,  0.974345}}}, {-1.62983, -1.26481, -2.54835, -1.71152, -1.47398,  -3.38347}}, {{{0.698099, 0.622589, 0.445884, 0.0971072, 0.0235963,  0.193392, 0.49464, 0.00703223, 0.453891, 0.185926, 0.922717,  0.151263, 0.401193, 0.0662702, 0.919217, 0.0569116, 0.179008,  0.196106, 0.563693, 0.21306, 0.27191, 0.436875, 0.631537,  0.238715, 0.573811, 0.814286, 0.185653, 0.141608, 0.550214,  0.620894, 0.691013, 0.134576, 0.0963236, 0.434969, 0.768296,  0.983313, 0.69513, 0.368699, 0.849078, 0.926401, 0.516123,  0.172593, 0.285385, 0.713341, 0.244213, 0.735718, 0.653848,  0.474625, 0.670402, 0.921431, 0.468195, 0.333017, 0.120188,  0.300537, 0.777183, 0.198441, 0.0238644, 0.865568, 0.0088872,  0.215129, 0.328734, 0.496869, 0.159809, 0.288728, 0.812611,  0.324277, 0.874424, 0.575387, 0.568398, 0.588559, 0.220575,  0.100761, 0.897996, 0.667128, 0.75238, 0.767744, 0.777808,  0.366591, 0.975198, 0.569303, 0.753943, 0.501023, 0.96631,  0.354174,  0.425209}, {{2, 0.00415373}, {4, 0.806501}, {0, 0.0654465}, {10,  0.612598}}}, {-1.79063, -1.74854, -1.73245, -1.96299, -1.40202,  -2.34626}}, {{{0.679877, 0.932078, 0.49006, 0.0441999, 0.0913179,  0.711502, 0.389299, 0.146204, 0.42419, 0.959122, 0.621555,  0.368396, 0.0575988, 0.983924, 0.052252, 0.614453, 0.556576,  0.0176142, 0.698078, 0.189243, 0.552422, 0.211113, 0.632631,  0.576645, 0.872545, 0.279035, 0.142571, 0.532445, 0.781227,  0.567533, 0.753273, 0.386241, 0.357037, 0.608411, 0.131718,  0.0178451, 0.299438, 0.624487, 0.079466, 0.403392, 0.742863,  0.606872, 0.381388, 0.214149, 0.190441, 0.395759, 0.748757,  0.637504, 0.317896, 0.116724, 0.606185, 0.105059, 0.536669,  0.549191, 0.852912, 0.718817, 0.179632, 0.94078, 0.721194,  0.700972, 0.880193, 0.316294, 0.641728, 0.29758, 0.137331,  0.709421, 0.26034, 0.0834308, 0.94689, 0.313662, 0.511584,  0.445927, 0.628995, 0.196938, 0.905398, 0.340868, 0.0923261,  0.647747, 0.0524861, 0.622051, 0.912695, 0.706966, 0.331292,  0.921079,  0.0325013}, {{2, 0.390673}, {12, 0.689563}, {9, 0.623499}, {0,  0.89517}}}, {-1.24188, -1.28377, -1.64583, -2.11965, -2.67738,  -2.94701}}, {{{0.681252, 0.429223, 0.540068, 0.94828, 0.36759,  0.917639, 0.094141, 0.319285, 0.170652, 0.0122407, 0.753273,  0.226959, 0.522905, 0.959755, 0.131222, 0.314265, 0.815939,  0.628463, 0.210143, 0.281764, 0.425266, 0.9389, 0.586644,  0.386593, 0.744015, 0.509677, 0.0465764, 0.438313, 0.376425,  0.592038, 0.952435, 0.119027, 0.205773, 0.579797, 0.199163,  0.892068, 0.682868, 0.620042, 0.0679411, 0.577803, 0.866929,  0.991579, 0.857798, 0.29604, 0.441663, 0.0526795, 0.271154,  0.909447, 0.697648, 0.543003, 0.224577, 0.471134, 0.321223,  0.950965, 0.272142, 0.352106, 0.11545, 0.371168, 0.0729792,  0.460038, 0.432583, 0.751126, 0.00503812, 0.882235, 0.565654,  0.759547, 0.14724, 0.586196, 0.123991, 0.706867, 0.876086,  0.676749, 0.426343, 0.163864, 0.651509, 0.205615, 0.10512,  0.212899, 0.379367, 0.853509, 0.98967, 0.841731, 0.306388,  0.39347,  0.557088}, {{7, 0.0906055}, {14, 0.301349}, {14, 0.511235}, {13,  0.991434}}}, {-1.18895, -1.97695, -1.75092, -2.47248, -1.45215,  -2.73472}}, {{{0.331059, 0.154109, 0.925039, 0.867442, 0.624192,  0.278023, 0.248291, 0.441099, 0.460327, 0.626514, 0.0426755,  0.335978, 0.247428, 0.247148, 0.189167, 0.346308, 0.405697,  0.94076, 0.795697, 0.789221, 0.315091, 0.639411, 0.284461,  0.797787, 0.984032, 0.485301, 0.359422, 0.930345, 0.35984,  0.207278, 0.111131, 0.489246, 0.899513, 0.580764, 0.0684559,  0.153268, 0.652085, 0.333616, 0.879289, 0.806959, 0.246388,  0.392856, 0.0835926, 0.0177388, 0.931297, 0.753445, 0.799131,  0.219952, 0.947265, 0.268144, 0.439709, 0.289607, 0.587425,  0.0608661, 0.328578, 0.800361, 0.687912, 0.480103, 0.260122,  0.647093, 0.0358271, 0.146487, 0.380833, 0.840134, 0.789439,  0.753631, 0.29724, 0.822395, 0.858141, 0.000185174, 0.498109,  0.602443, 0.910876, 0.732041, 0.0583998, 0.312836, 0.323451,  0.671175, 0.729822, 0.512476, 0.635539, 0.191072, 0.4697,  0.865383,  0.599712}, {{9, 0.0445858}, {15, 0.0888673}, {1, 0.025249}, {6,  0.810273}}}, {-2.03362, -1.72501, -1.38264, -1.60279, -2.0633,  -2.19208}}, {{{0.290955, 0.791627, 0.202854, 0.952132, 0.29077,  0.293518, 0.600411, 0.0412557, 0.558729, 0.235118, 0.287575,  0.717805, 0.887554, 0.505296, 0.775099, 0.0822655, 0.696482,  0.0355962, 0.909716, 0.482554, 0.651896, 0.946729, 0.884467,  0.67228, 0.360941, 0.155102, 0.681613, 0.720149, 0.0701707,  0.861584, 0.0812024, 0.678893, 0.511442, 0.626465, 0.793628,  0.961088, 0.623887, 0.121169, 0.0185287, 0.878823, 0.927405,  0.0855728, 0.108812, 0.396269, 0.275509, 0.138844, 0.224345,  0.723989, 0.914569, 0.983742, 0.542731, 0.00384031, 0.844398,  0.122158, 0.461529, 0.324947, 0.332956, 0.495693, 0.667901,  0.363859, 0.709069, 0.374524, 0.649373, 0.485036, 0.781664,  0.288951, 0.54056, 0.0887667, 0.506154, 0.150107, 0.316216,  0.364778, 0.591586, 0.166365, 0.773484, 0.360937, 0.747188,  0.0442071, 0.311955, 0.0359901, 0.414232, 0.548514, 0.644054,  0.672131,  0.705162}, {{14, 0.17399}, {15, 0.994681}, {11, 0.187095}, {9,  0.923499}}}, {-2.11107, -1.63796, -1.80587, -1.40532, -1.56718,  -2.71412}}, {{{0.885039, 0.454121, 0.0983283, 0.417344, 0.734931,  0.137905, 0.733551, 0.825759, 0.568566, 0.364421, 0.372613,  0.0785709, 0.524359, 0.0524658, 0.336623, 0.664339, 0.975845,  0.408412, 0.664492, 0.959177, 0.801855, 0.413731, 0.477397,  0.035678, 0.916816, 0.95961, 0.379068, 0.618334, 0.181884,  0.821705, 0.645518, 0.792575, 0.613318, 0.457284, 0.272905,  0.714004, 0.0889594, 0.404818, 0.936282, 0.0496646, 0.113115,  0.996406, 0.27179, 0.0904878, 0.31126, 0.582675, 0.794393,  0.0548097, 0.394444, 0.623065, 0.415325, 0.436476, 0.21256,  0.801361, 0.769807, 0.643901, 0.599242, 0.344077, 0.496902,  0.929897, 0.510282, 0.939258, 0.560621, 0.880233, 0.397168,  0.942852, 0.28883, 0.789745, 0.0859076, 0.360177, 0.494437,  0.734935, 0.691463, 0.737111, 0.0791118, 0.298459, 0.478903,  0.935751, 0.309305, 0.654558, 0.879662, 0.591674, 0.812402,  0.72466,  0.369379}, {{8, 0.652416}, {12, 0.251782}, {0, 0.844428}, {15,  0.972212}}}, {-1.81144, -1.3156, -1.73998, -2.11744, -1.49119,  -3.05123}}, {{{0.709564, 0.962951, 0.0546827, 0.886304, 0.349387,  0.468514, 0.319747, 0.194841, 0.612276, 0.389402, 0.0212883,  0.715937, 0.676526, 0.0800978, 0.36673, 0.836276, 0.0848515,  0.267696, 0.64207, 0.466896, 0.432436, 0.0159142, 0.797642,  0.494685, 0.722872, 0.0529631, 0.74296, 0.60838, 0.373485,  0.584449, 0.423212, 0.41354, 0.761208, 0.195047, 0.401924,  0.697602, 0.0846829, 0.114949, 0.0351937, 0.861327, 0.999831,  0.847253, 0.393124, 0.39443, 0.567396, 0.831339, 0.595481,  0.899746, 0.844524, 0.778376, 0.852521, 0.291365, 0.471039,  0.193927, 0.429309, 0.877826, 0.70983, 0.99888, 0.0273847,  0.180223, 0.625148, 0.883932, 0.992191, 0.318896, 0.625316,  0.0366785, 0.599067, 0.924466, 0.0579208, 0.20534, 0.00358633,  0.0247202, 0.213397, 0.426964, 0.151065, 0.733355, 0.742358,  0.233037, 0.721756, 0.855529, 0.0325278, 0.234156, 0.694371,  0.675306,  0.40738}, {{8, 0.350225}, {2, 0.70218}, {9, 0.35641}, {4,  0.782064}}}, {-1.8528, -1.82827, -1.29519, -1.97015, -1.55342,  -2.8524}}, {{{0.313546, 0.103113, 0.431944, 0.724143, 0.108207,  0.0995265, 0.407223, 0.510746, 0.681243, 0.948462, 0.673868,  0.768388, 0.448206, 0.226705, 0.818339, 0.73586, 0.21405,  0.532334, 0.143033, 0.32848, 0.863825, 0.830154, 0.786623,  0.546416, 0.550278, 0.727041, 0.35468, 0.822273, 0.442072,  0.627514, 0.947456, 0.311527, 0.760829, 0.679053, 0.273588,  0.543139, 0.312623, 0.452347, 0.455249, 0.807279, 0.098573,  0.920013, 0.312216, 0.4788, 0.234748, 0.0898596, 0.525592,  0.932384, 0.68447, 0.362819, 0.170912, 0.110111, 0.242398,  0.735305, 0.223456, 0.798584, 0.48157, 0.0562518, 0.949868,  0.255444, 0.168947, 0.603904, 0.494619, 0.448165, 0.070374,  0.683891, 0.182404, 0.969365, 0.835626, 0.594032, 0.656811,  0.0369817, 0.151156, 0.231213, 0.485899, 0.926871, 0.908758,  0.495908, 0.262443, 0.128287, 0.427188, 0.439657, 0.312575,  0.872843,  0.258241}, {{16, 0.835752}, {0, 0.817955}, {5, 0.424678}, {5,  0.187867}}}, {-1.70133, -2.09112, -1.64194, -1.16758, -2.07145,  -2.76016}}}; 

    for (unsigned int i = 0; i < sizeof (data) / sizeof (data[0]); ++i)
      {
        LogisticRegression r (17, 6, 1, data[i].inputs.w, 1);
        std::vector<Feature> features (data[i].inputs.f, data[i].inputs.f + 4);
        IndexedFeatureSet feature_set;
        unsigned int ns = 256 * drand48 ();
        feature_set[ns] = features;

        std::vector<float> logp = r.logp (feature_set);
  
        for (unsigned int j = 0; j < 5; ++j)
          {
            assert (fabs (logp[j] - data[i].logp[j]) <=
                    1e-4 * (1 + fabs (logp[j]) + fabs (data[i].logp[j])) ||
                    (std::cerr << "logp[" << j << "] " << logp[j] 
                               << " ?= " << data[i].logp[j] << std::endl,
                     0));
          }
      }
  }

  struct GradWTestData
  {
    struct
      {
        float w[9 * 3];
        Feature f[2];
        float   q[4];
      }                 inputs;
  
    struct
      {
        float logp[4];
        float kl;
        float gradw[9 * 3];
      }                 outputs;
  };

//  std::ostream&
//  operator<< (std::ostream& o,
//              const std::vector<float>& v)
//    {
//      for (unsigned int i = 0; i < v.size (); ++i)
//        {
//          o << v[i] << ", ";
//        }
//
//      return o;
//    }

  void
  test_gradw (void)
  {
    GradWTestData data[] = 
    {{{{0.905638, 0.442821, 0.872684, 0.558764, 0.890784, 0.380951,  0.31314, 0.842586, 0.565679, 0.95908, 0.670511, 0.579514,  0.541265, 0.822774, 0.434008, 0.154929, 0.9386, 0.612426,  0.332957, 0.517111, 0.580024, 0.049592, 0.532275, 0.290892,  0.674385, 0.606771,  0.659591}, {{2, 0.732128}, {7, 0.783601}}, {0.134438, 0.206254,  0.529572,  0.129736}}, {{-1.68067, -0.914815, -1.30669, -1.94879}, -1.35945,  {0, 0, 0, 0, 0, 0, -0.0379327, -0.14228, 0.189518, 0, 0, 0, 0, 0, 0,  0, 0, 0, 0, 0, 0, -0.0405996, -0.152283, 0.202842, 0, 0,  0}}}, {{{0.26674, 0.675941, 0.310028, 0.676657, 0.443966,  0.241933, 0.155099, 0.738057, 0.83154, 0.908976, 0.637988,  0.158033, 0.781948, 0.376702, 0.347096, 0.483648, 0.175177,  0.71711, 0.614968, 0.700047, 0.949357, 0.370659, 0.725426,  0.482124, 0.682617, 0.694718,  0.415398}, {{6, 0.805467}, {2, 0.238651}}, {0.381259, 0.21918,  0.0567612,  0.3428}}, {{-1.47166, -1.26401, -1.04089, -2.00401}, -1.58418, {0,  0, 0, 0, 0, 0, 0.0362067, -0.0151161, -0.0707314, 0, 0, 0, 0, 0,  0, 0, 0, 0, 0.122201, -0.0510181, -0.238725, 0, 0, 0, 0, 0,  0}}}, {{{0.543809, 0.622312, 0.909377, 0.625163, 0.167107,  0.275215, 0.425729, 0.449986, 0.449997, 0.660247, 0.725682,  0.500629, 0.0793378, 0.934821, 0.243558, 0.818013, 0.38462,  0.519423, 0.438091, 0.579362, 0.931835, 0.259123, 0.370681,  0.172251, 0.388026, 0.636812,  0.461304}, {{1, 0.547088}, {3, 0.220919}}, {0.285802, 0.0281187,  0.0767488,  0.60933}}, {{-1.16338, -1.39952, -1.3901, -1.65126}, -1.4847, {0,  0, 0, -0.0145662, -0.119591, -0.0942645, 0, 0,  0, -0.00588198, -0.048292, -0.0380649, 0, 0, 0, 0, 0, 0, 0, 0, 0,  0, 0, 0, 0, 0, 0}}}, {{{0.70135, 0.309894, 0.596473, 0.691585,  0.766529, 0.0663363, 0.77846, 0.306965, 0.247106, 0.628245,  0.199098, 0.37513, 0.987983, 0.257565, 0.0268468, 0.987103,  0.351171, 0.79626, 0.479758, 0.766184, 0.989574, 0.760685,  0.382656, 0.995261, 0.288225, 0.450791, 0.786183}, {{6, 0.303676},  {2, 0.521696}}, {0.231115, 0.00464238, 0.599173,  0.16507}}, {{-1.19807, -1.35707, -1.32046, -1.74988}, -1.36322,  {0, 0, 0, 0, 0, 0, -0.0368639, -0.131871, 0.173286, 0, 0, 0, 0, 0, 0,  0, 0, 0, -0.0214583, -0.0767611, 0.100869, 0, 0, 0, 0, 0,  0}}}, {{{0.756209, 0.808624, 0.621582, 0.286607, 0.498644,  0.781778, 0.634478, 0.935436, 0.702384, 0.302019, 0.868294,  0.945862, 0.9417, 0.919363, 0.873033, 0.657637, 0.490909, 0.13318,  0.569356, 0.135941, 0.106455, 0.125458, 0.572645, 0.861351,  0.350246, 0.316833,  0.951063}, {{4, 0.574744}, {6, 0.851601}}, {0.326221, 0.19302,  0.389782,  0.0909768}}, {{-0.989194, -1.37113, -1.42287, -2.01529}, -1.3253,  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.0262403, -0.0349446,  0.0854991, 0, 0, 0, -0.0388804, -0.0517776, 0.126684, 0, 0, 0, 0,  0, 0}}}, {{{0.233037, 0.448291, 0.693445, 0.207517, 0.313674,  0.575258, 0.0358082, 0.716608, 0.180493, 0.00590204, 0.899867,  0.610154, 0.0550357, 0.433257, 0.0385158, 0.259908, 0.738202,  0.482194, 0.463772, 0.408307, 0.203146, 0.16561, 0.824465,  0.25909, 0.97011, 0.717319,  0.131019}, {{6, 0.0515726}, {8, 0.656436}}, {0.135531, 0.0908347,  0.319568,  0.454066}}, {{-1.07498, -1.24378, -1.63923, -1.73571}, -1.57065,  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,  0, -0.0106123, -0.0101834, 0.00646921, 0, 0,  0, -0.135077, -0.129618, 0.0823426}}}, {{{0.136158, 0.195344,  0.724811, 0.420907, 0.702901, 0.156828, 0.464902, 0.682705,  0.220707, 0.693056, 0.0565955, 0.479558, 0.0550971, 0.868592,  0.797506, 0.509448, 0.337778, 0.737573, 0.745933, 0.853012,  0.195718, 0.642362, 0.410969, 0.37707, 0.0595596, 0.447018,  0.686158}, {{1, 0.956163}, {7, 0.356659}}, {0.31513, 0.240272,  0.296961,  0.147637}}, {{-1.2367, -1.0496, -1.58382, -1.86826}, -1.38807, {0,  0, 0, 0.0237033, -0.104993, 0.0877491, 0, 0, 0, 0, 0, 0, 0, 0, 0,  0, 0, 0, 0, 0, 0, 0.00884157, -0.0391635, 0.0327313, 0, 0,  0}}}, {{{0.597133, 0.16466, 0.7939, 0.0808549, 0.728542, 0.367155,  0.284452, 0.743076, 0.990969, 0.621222, 0.431439, 0.547358,  0.348608, 0.210253, 0.0543696, 0.487799, 0.90159, 0.524095,  0.0982065, 0.13114, 0.6114, 0.302839, 0.824748, 0.995188,  0.0142669, 0.138179,  0.0308476}, {{5, 0.914333}, {6, 0.285725}}, {0.388732, 0.376315,  0.0863437,  0.148609}}, {{-1.45702, -1.06926, -1.2772, -1.93109}, -1.36602,  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.142455,  0.0302221, -0.175985, 0.0445165, 0.00944426, -0.0549947, 0, 0, 0,  0, 0, 0}}}, {{{0.149803, 0.314956, 0.623899, 0.946148, 0.93955,  0.260587, 0.1361, 0.0445583, 0.415455, 0.16238, 0.00495959,  0.433158, 0.112616, 0.337632, 0.00977139, 0.418891, 0.974437,  0.306784, 0.0954381, 0.133166, 0.203413, 0.560389, 0.924181,  0.83841, 0.0536099, 0.245433,  0.300283}, {{6, 0.892262}, {6, 0.11406}}, {0.365389, 0.0609136,  0.314507,  0.25919}}, {{-1.40165, -1.36369, -1.29299, -1.49769}, -1.39006,  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,  0.119952, -0.196035, 0.0403128, 0, 0, 0, 0, 0, 0}}}, {{{0.822466,  0.159223, 0.414546, 0.585989, 0.484834, 0.149452, 0.995655,  0.611552, 0.17805, 0.0540137, 0.862489, 0.408139, 0.617661,  0.129833, 0.0240787, 0.354529, 0.372228, 0.82955, 0.131817,  0.240469, 0.387382, 0.665367, 0.284113, 0.541864, 0.564916,  0.506144, 0.869568}, {{2, 0.955876}, {7, 0.0800816}}, {0.144004,  0.352817, 0.139011,  0.364169}}, {{-0.912731, -1.31042, -1.70415, -1.91774}, -1.52905,  {0, 0, 0, 0, 0, 0, -0.246064, 0.0794419, -0.0410226, 0, 0, 0, 0, 0, 0,  0, 0, 0, 0, 0, 0, -0.0206148, 0.0066555, -0.0034368, 0, 0, 0}}}}; 

    for (unsigned int i = 0; i < sizeof (data) / sizeof (data[0]); ++i)
      {
        LogisticRegression rzero (9, 4, 1, data[i].inputs.w, 1);
        std::vector<Feature> features (data[i].inputs.f, data[i].inputs.f + 2);
        IndexedFeatureSet feature_set;
        unsigned int ns = 256 * drand48 ();
        feature_set[ns] = features;
        std::vector<float> logp = rzero.logp (feature_set);
        std::vector<float> q (data[i].inputs.q, data[i].inputs.q + 4);
        float kl = kldivergence (logp, q);

        assert (fabs (kl - data[i].outputs.kl) <=
                1e-4 * (1 + fabs (kl) + fabs (data[i].outputs.kl)) ||
                (std::cerr << "kl = " << kl
                           << " ?= " << data[i].outputs.kl << std::endl,
                 0));

        for (unsigned int j = 0; j < 4; ++j)
          {
            assert (fabs (logp[j] - data[i].outputs.logp[j]) <=
                    1e-4 * (1 + fabs (logp[j]) + fabs (data[i].outputs.logp[j])) ||
                    (std::cerr << "logp[" << j << "] " << logp[j] 
                               << " ?= " << data[i].outputs.logp[j] << std::endl,
                     0));
          }

        float xdotx = 0;

        for (unsigned int k = 0; k < 2; ++k)
          { 
            xdotx += features[k].value;
          }

        LogisticRegression::DFunc dfunc (logp, q, xdotx);
          
        PianoVec s0 (3, 0);
        PianoVec sprime0 (dfunc (0, s0));

        for (unsigned int k = 0; k < 3; ++k)
          {
            std::vector<float> dw (9, 0);

            for (unsigned int j = 0; j < 2; ++j)
              {
                dw[features[j].feature] += sprime0[k] * features[j].value;
              }

            for (unsigned int j = 0; j < 9; ++j)
              {
                assert (fabs (dw[j] - data[i].outputs.gradw[3 * j + k]) <=
                        1e-5 * (1 + fabs (dw[j]) + fabs (data[i].outputs.gradw[3 * j + k])) ||
                        (std::cerr << "(" << k << ") dw[" << j << "] " << dw[j] 
                                   << " ?= " << data[i].outputs.gradw[3 * j + k] << std::endl,
                         0));
              }
          }
      }
  }

  struct UpdateTestData
  {
    struct
      {
        float w[9 * 3];
        Feature f[2];
        float   q[4];
      }                 inputs;
  
    struct
      {
        float wnew[9 * 3];
      }                 outputs;
  };

  void
  test_update (void)
  {
    UpdateTestData data[] = 
{{{{0.0106965, 0.527527, 0.256434, 0.704891, 0.157111, 0.0785289, 
    0.23372, 0.156676, 0.770983, 0.512371, 0.560752, 0.207422, 
    0.549072, 0.900001, 0.0269025, 0.496306, 0.787705, 0.712925, 
    0.191926, 0.390111, 0.0746026, 0.609365, 0.0193986, 0.548841, 
    0.0639062, 0.0818387, 
    0.762964}, {{3, 0.843949}, {4, 0.906795}}, {0.00244149, 0.390402, 
    0.506973, 0.100183}}, {{0.0106965, 0.527527, 0.256434, 0.704891, 
    0.157111, 0.0785289, 0.23372, 0.156676, 0.770983, 0.300177, 
    0.52752, 0.471439, 0.321076, 0.864294, 0.310579, 0.496306, 
    0.787705, 0.712925, 0.191926, 0.390111, 0.0746026, 0.609365, 
    0.0193986, 0.548841, 0.0639062, 0.0818387, 
    0.762964}}}, {{{0.490938, 0.968492, 0.479851, 0.58674, 0.590938, 
    0.94159, 0.983544, 0.799034, 0.878013, 0.749664, 0.593433, 
    0.724432, 0.268647, 0.730265, 0.0445928, 0.660525, 0.186808, 
    0.967301, 0.200644, 0.75373, 0.183499, 0.438057, 0.513371, 
    0.617919, 0.69256, 0.469564, 
    0.03352}, {{7, 0.0311789}, {6, 0.101623}}, {0.510759, 0.0483461, 
    0.224575, 0.216319}}, {{0.490938, 0.968492, 0.479851, 0.58674, 
    0.590938, 0.94159, 0.983544, 0.799034, 0.878013, 0.749664, 
    0.593433, 0.724432, 0.268647, 0.730265, 0.0445928, 0.660525, 
    0.186808, 0.967301, 0.227299, 0.731944, 0.181015, 0.446235, 
    0.506686, 0.617156, 0.69256, 0.469564, 0.03352}}}, {{{0.778311, 
    0.456542, 0.507713, 0.954963, 0.0480458, 0.411949, 0.847188, 
    0.768154, 0.0807449, 0.211306, 0.0934574, 0.584656, 0.642688, 
    0.697935, 0.475539, 0.892095, 0.173124, 0.664415, 0.44436, 
    0.790473, 0.645149, 0.61444, 0.212215, 0.566863, 0.866838, 
    0.157897, 0.704502}, {{8, 0.6119}, {7, 0.818792}}, {0.234202, 
    0.269168, 0.264908, 0.231722}}, {{0.778311, 0.456542, 0.507713, 
    0.954963, 0.0480458, 0.411949, 0.847188, 0.768154, 0.0807449, 
    0.211306, 0.0934574, 0.584656, 0.642688, 0.697935, 0.475539, 
    0.892095, 0.173124, 0.664415, 0.44436, 0.790473, 0.645149, 
    0.514259, 0.280296, 0.520199, 0.791972, 0.208776, 0.66963}}}};

    for (unsigned int i = 0; i < sizeof (data) / sizeof (data[0]); ++i)
      {
        std::vector<float> w (data[i].inputs.w, data[i].inputs.w + 9 * 3);
        LogisticRegression rzero (9, 4, 1, &w[0], 1);
        std::vector<Feature> features (data[i].inputs.f, data[i].inputs.f + 2);
        IndexedFeatureSet feature_set;
        unsigned int ns = 256 * drand48 ();
        feature_set[ns] = features;
        std::vector<float> logp = rzero.logp (feature_set);
        std::vector<float> q (data[i].inputs.q, data[i].inputs.q + 4);

        rzero.update (1, 1, 0, feature_set, logp, q);
  
        assert (rzero.failures () == 0);

        for (unsigned int j = 0; j < 9 * 3; ++j)
          {
            assert (fabs (w[j] - data[i].outputs.wnew[j]) <=
                    1e-3 * (1 + fabs (w[j]) + fabs (data[i].outputs.wnew[j])) ||
                    (std::cerr << "(" << i << ") w[" << j << "] " << w[j] 
                               << " ?= " << data[i].outputs.wnew[j] << std::endl,
                     0));
          }
      }
  }
}

int
main (void)
{
  test_logp ();
  test_gradw ();
  test_update ();

  return 0;
}
